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We propose a two-dimensional model for a multilayer stack of dipolar Bose-Einstein condensates 
formed by a strong optical lattice. We derive effective intra- and interlayer dipole-dipole interac- 
tion potentials and provide simple analytical approximations for a given number of lattice sites at 
arbitrary polarization. We find that the interlayer dipole-dipole interaction changes the transverse 
aspect ratio of the ground state in the central layers depending on its polarization and the number 
of lattice sites. The changing aspect ratio should be observable in time of flight images. Further- 
more, we show that the interlayer dipole-dipole interaction reduces the excitation energy of local 
perturbations affecting the development of a roton minimum. 
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I. INTRODUCTION 

Layered structures of magnetic materials play a cru- 
cial role both in today's technology and in fundamen- 
tal physical theories. Technological examples are aplenty 
in the magneto-electronic industries, e.g., hard disks or 
magnetic sensors. One theoretical goal of studying mul- 
tilayers is to illuminate the elusive theory of high-T c su- 
perconductivity, where the layered structure appears to 
play a crucial role [1]. For a realistic theory of atomic or 
molecular multilayers it is, however, vital to include the 
dipole-dipole interaction (DDI) between the underlying 
particles. 

The study of magnetic single- and multilayer films has 
enjoyed a long history in condensed matter physics (for a 
recent review, see Ref. [2] and references therein). There, 
an alternating structure of ferromagnetic and nonmag- 
netic layers is deposited on a substrate, e.g., by atomic 
beam epitaxy. However, structural instabilities induced, 
e.g., by temperature changes and film thickness variation 
often complicate experiments in thin films. 

Quantum-degenerate dipolar gases have received much 
attention recently from both theoretical and experimen- 
tal studies (for recent reviews, see Refs. [3, 4]). Their 
DDI crucially affects the ground-state properties [5, 6], 
stability [7-9], and dynamics of the gas [10]. Further- 
more, they offer a route for studying exciting many- 
body quantum effects, such as a superfluid-to-crystal 
quantum phase transition [11], supersolids [12] or even 
topological order [13]. Recent advances in experimen- 
tal techniques have paved the way for a Bose-Einstein 
condensate (BEC) of 52 Cr with a magnetic dipole mo- 
ment 6fiB (Bohr magneton hb), much larger than con- 
ventional alkali BECs [14-16]. Promising candidates for 
future dipolar BEC experiments are Er and Dy with 
even larger magnetic moments of 1[Lb and 10/is, respec- 
tively [17, 18]. Furthermore, DDI-induced decoherence 
and spin textures have been observed in alkali-metal con- 
densates [19, 20]. Dipolar effects also play a crucial role 
in experiments with Rydberg atoms [21] and heteronu- 
clear molecules [22, 23]. Bosonic heteronuclear molecules 
may provide a basis for future experiments on BECs with 



dipole moments much larger than in atomic BECs [24] . 

In contrast to solid state thin film structures, the layer 
width and spacing of BECs in optical lattices are pre- 
cisely tunable with external fields. This makes dipolar 
BECs a prime candidate for investigating the effects of 
DDI in multilayers. For example, it has been shown that 
the DDI stabilizes quasi-two-dimensional ultracold gases 
for perpendicular polarization [9, 25] and enables con- 
trolled chemical reactions [23]. Another intriguing effect 
is the occurrence of interlayer bound states [26-31]. How- 
ever, it is still unclear to what extend effective models for 
multilayers of dipolar BEC at arbitrary polarization are 
valid and how interlayer DDI can be detected. 

In this article, we investigate the effect of interlayer 
DDI on the ground state of the BEC. We present an 
effective two-dimensional (2D) model for an arbitrarily 
polarized dipolar BEC in a strong one-dimensional (ID) 
optical lattice. Our 2D model offers a clear advantage for 
numerical computation of ground state properties com- 
pared to computations for a full three-dimensional (3D) 
Gross-Pitaevskii equation (GPE) : our computation times 
reduce to seconds instead of dozens of hours. Previously, 
such dimension-reduced models have been derived for 
BECs without DDI [32-38] and with dipolar interactions 



x 




FIG. 1. (Color online) Setup of the multilayered dipolar BEC 
polarized along d. An optical lattice along z separates the 
dipolar BEC into 2D layers in the x-y plane with distance 5. 
Apart from the intralayer DDI U2U, each layer interacts with 
other layers via the interlayer DDI f/| D . 
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in a single layer [39, 40]. We also derive the effective 2D 
intra- and interlayer DDI potentials governing the layers 
of quasi-2D BECs. These potentials allow for useful an- 
alytical approximations, which were used in a previous 
work on multilayer dipolar BECs with perpendicular po- 
larization [27]. We establish that the 2D model is valid 
by comparing its ground states to ground states of the 
3D GPE for weakly interacting BECs at zero tempera- 
ture [41]. We suggest that the interlayer DDI is observ- 
able in the transverse aspect ratio of the central layers 
after time of flight expansion. Moreover, we calculate 
the Bogoliubov excitation energies for a transversely ho- 
mogeneous BEC with contact, intra- and interlayer DDI. 
The interlayer DDI reduces the squared Bogoliubov en- 
ergy and, therefore, influences the occurance of a roton 
minimum. 

In Sec. II we present our 2D model and effective intra- 
and interlayer potentials for a dipolar BEC trapped in a 
strong ID optical lattice. We also present a single mode 
approximation valid for the central layers of the BEC. 
In Sec. Ill we compare ground states of our model and 
its single mode approximation to ground states of the 
3D GPE. We find good agreement between these ground 
states, which indicates the validity of our model. In 
Sec. IV we compute numerically the aspect ratio of the 
BEC in the central layer as a function of the number of 
lattice sites and polarization direction. We find a marked 
change in the aspect ratio owing to the interlayer DDI, 
which should be observable in experiments. In Sec. V 
we derive the Bogoliubov dispersion for a transverse ho- 
mogeneous, multilayered dipolar BEC. We conclude in 
Sec. VI. In App. A we give a detailed derivation of the 
2D model presented in Sec. II. 



II. EFFECTIVE 2D MODEL 

We consider a dilute dipolar BEC at zero temperature 
trapped in a transverse harmonic potential Vii (x, y) — 
^^^-(x 2 + y 2 ) and a longitudinal optical lattice V (z) = 
Vq sin 2 (/c;z). Here, m is the particle mass, lu the trap 
frequency, Vq the lattice height, and ki the wave num- 
ber of the lattice laser. We focus on atomic BECs 
with a magnetic dipole moment but it is straightfor- 
ward to extend the analysis to degenerate bosonic gases 
with electric dipole moments. We assume that an ex- 
ternal field polarizes the atoms along a normalized axis 
d = (d x ,d y ,d z ) = (cos (fismd, sin^>sin$, cosd) with (f> 
and $ the azimuthal and polar angles, respectively. Then 
the dipole-dipole interaction (DDI) is described by 
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where Cdd = Mo-D 2 with (iq is the magnetic vacuum per- 
meability and D the dipole moment (for electric dipoles 
Cdd = D 2 /eo, where eo is the vacuum permittivity). We 
note that it is possible to modify the DDI strength c<jd 
by means of a rotating magnetic field [42]. 



At zero temperature, a weakly interacting BEC is de- 
scribed by the GPE [41]. For simplicity, we introduce di- 
mensionless quantities by rescaling lengths with the lat- 
tice distance 8 = ir/ki, that is, r — > r<5, energies with 
H 2 /mS 2 = 2E r /ir 2 (E r is the recoil energy), and the 
wave function of the gas with the central density n(0), 
if) — > i/jy/n(0). In these units the normalization of the 
wave function is J <i 3 r|V>(r, t)\ 2 = N/n(0)S 3 with N the 
total number of atoms. Away from shape resonances, the 
wave function i/j = ip(r, t) of the dipolar BEC is governed 
by the GPE [6, 43, 44] 
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Here, g — 4ira s n(0)5 2 is the dimensionless contact inter- 
action strength with a s the s-wave scattering length and 
9d = mc dd n(Q)8 2 /Zh 2 is the dimensionless DDI strength. 
Furthermore, V\ lo {p) = (m 2 w 2 5 A /2H 2 )p 2 with p — (x,y) 
and V (z) — {Vqtt 2 /2) sin 2 (7rz), where Vq is the lattice 
amplitude in units of the recoil energy E r . The nonlocal 
dipolar potential V^d is given by 

V dd (r) = -3g d d dd J d 3 r'U 3D (r- r')|V(r',*)| 2 (3) 

with the kernel Usi)(r) = l/47r|r| and the notation d d = 
d • V, d dd = d 2 



A. Coupled modes 

For strong optical lattices we derive an effective 2D 
equation for the wave function on each lattice site. This 
is possible because a strong optical lattice with V) hu 
causes the BEC to form layers separated by the lat- 
tice distance 6 (cf. Fig. 1) [41, 45]. We assume that 
the axial extend 7 of the BEC in each layer is much 
larger than the s-wave scattering length. Additionally, 
in the quasi-2D regime 7~ 2 3> \g — g d \ [46]. This con- 
dition allows us to approximate the optical lattice as a 
train of harmonic potentials and the axial wave func- 
tion as its ground state. Then the wave function sepa- 
rates into ip(r,t) = e- lt / 2 ^ 2 Y,ti>t(p,t)m{z) [33, 40, 41]. 
The sum extends over all lattice sites I. Under our as- 
sumptions the axial wave function on each site i at posi- 
tion is described by a Gaussian wi(z) — w(z — zi) = 
(l/7r7 2 ) 1 / 4 e _ ( z_z ^ I 21 ; the Gaussians do not mutually 
overlap (J dzwe(z)wj(z) ~ for £ =^ j). In the quasi- 
2D limit 7~ 2 = \/Vqit 2 . More generally, in a homoge- 
neous BEC it is also possible to treat the layer width 
7 as a variational parameter that minimizes the Gross- 
Pitaevskii energy functional [9]. By inserting this wave 
function into Eq. (2) and integrating out the z direction 
we obtain the following equation for the radial wave func- 
tion tp£ = ipe{p, t) at site £ 
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FIG. 2. (Color online) Cumulative interlayer DDI C/| D for 
52 Cr at different layer separations \£ — j\. The solid lines show 
the interlayer DDI [Eq. (9)], whereas the dashed line shows 
the approximation Eq. (10) for nearest neighbors [Eq. (10) 
is indistinguishable from the solid lines for larger distances]. 
The dotted line indicates the intralayer DDI. The inset shows 
the intralayer DDI and the rectangle within indicates the ex- 
tend of the main panel. We set Vo = 30£V- 



Here, g — g/\/2ir"f and cjd — gd/\2nj are the effective 
2D interactions strengths. In the remainder of this article 
we neglect strongly suppressed terms in the effective DDI 
potential V^ D (see Appendix A for details). We find the 
following expression for its Fourier transform V^ D (k) = 
J^V^rjKk) with k = /c(cos</?, sin^) 

Kf D (k) = 3g d ]T ( [(4 cos <p + d y sin ^) 2 - d 2 } U^k) 

3 

+ 2id z (d x cos ip + d y sin tp)U^ dd (k)j \tpj\ 2 (k). 



(5) 



Here, 



" 2 7 2 



U: 



odd(fc) = ^ ^ 



7 2 fc + S e 
V2 7 



7 2 fc + S e 



j 2 k 



V2j 



-i 2 k - S, 



V2j 



(6) 



(7) 



where dij = (£ — j), i](x) = exp(a; 2 ) erfc(ir) and erfc(ir) = 
1 — erf (x) is the complementary error function. 

The effective dipolar interaction V^ D [Eq. (5)] contains 
both an intralayer DDI and an interlayer DDI. The in- 
tralayer DDI are the terms in Eq. (5) with I = j. By 
setting £ — j in Eqs. (6)-(7) we find that each layer expe- 
riences the effective DDI potential of a quasi-2D dipolar 
BEC [40]. The interlayer DDI are the terms in Eq. (5) 
with £ ^ j. For perpendicular polarization (d z = 1, 
d x = d y — 0) we recover the interlayer DDI potential 
discussed, e.g., in Ref. [27]. If the layer distance is much 



larger than the layer width (|<%| >• 7), i](x —> +00) van- 
ishes and the moduli of the kernels |f7evonl an d l^oddl be- 
come identical. Because of our assumption that 5 3> 7, 
this is fulfilled for the interlayer DDI between any two 
distinct sites. As a consequence, we split the total effec- 
tive DDI potential into a sum of intralayer and interlayer 
terms 

Vf u (k) = 3g d [(d x cos tp + dy simp) 2 - ^]f7 2D (fc)f^(k) 
+ 3g d ^2\dx cos ip + d y sin ip - id z sgn(%)] 2 

xt/^(fc)|^(k), 

(8) 



where sgn(x) is the sign of x. The kernels of this potential 
are Z7 2 d 



2C/ 2 °o and 
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In the limit of negligible layer width (7 <§; |<%|) the in- 
terlayer DDI in Eq. (9) can be approximated by 



-\Stj\k 



(10) 



This approximation becomes an identity in the limit 
7 — s- and nonzero \Sej\. The second line of Eq. (8) 
is the interlayer DDI potential for arbitrary polarization 
direction. Inserting approximation (10) into Eq. (8) for 
perpendicular polarization, we recover the interlayer DDI 
potential used in Refs. [27, 29]. We expect our general- 
ized interlayer DDI potential to be valid for bosons as 
well as fermions because fcrmions in different layers oc- 
cupy different quantum states. 

The kernel of the interlayer DDI potential IJ^(k) is 
shown in Fig. 2 as a cumulative plot over the five nearest 
lattice sites. For comparison we also show the intralayer 
DDI. Although not shown in Fig. 2, we established that 
for realistic parameters the potentials Ui^cn an d ^odd (^ or 
£ 7^ j) are indistinguishable from U^ D at the plot reso- 
lution. For interlayer interactions beyond nearest neigh- 
bors the approximation for f)| D in Eq. (10) becomes in- 
distinguishable from Eq. (9) . The interlayer DDI is linear 
in momentum for long wavelengths and drops exponen- 
tially for short wavelengths. It has been shown that this 
behavior leads to very weakly bound states in bilayer 
systems [27, 30, 31, 47]. According to Eq. (8) its sign is 
determined by the polarization direction. The interlayer 
and intralayer DDI for predominantly perpendicular po- 
larization < 7r/4) is attractive in momentum space for 
all k, whereas the interlayer DDI for predominantly par- 
allel polarization ($ > 7r/4) becomes repulsive for some 
k around the major axis with ip = <p. 
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B. Single mode approximation 

If we assume that the the BEC densities in each layer 
vary little over the central sites, we can simplify the 2D 
model to a single equation for the central site wave func- 
tion ipo(p). This assumption is reasonable for large lat- 
tices and we will test its validity in Sec. III. The single 
wave function ipo(p) approximates the wave functions in 
all lattice sites far from the boundaries. Consequently, we 
replace the effective dipolar potential V2 D (k) [Eq. (8)] by 
the site-local potential 

V 2 u(X) = 3g d ([(d x cos ip + dysmtp) 2 - d 2 z ]U 2 v(k) 

+ ^^[d x cos<p + d y s'mip — \d z sg\i(j)] 2 ij^{k) 



x |^ | 2 (k). 



(11) 



Inserting the inverse Fourier transform of Eq. (11) into 
Eq. (4) we are left with the uncoupled equation 



-2 V 2 + V ho + [g - g d (l - 3d 2 )] |0 O | 2 + U 2D 



(12) 

for the central site wave function ipo = ipo(p)- We assume 
a lattice that is symmetric around the central site so that 
the dipole terms linear in d z in Eq. (11) vanish after 
summation. Using Eq. (10) for U^n we can perform the 
summation in Eq. (11) and find 



V 2 D(k) ~ 3g d [(d x costp + d y simp) 2 - d 2 z 



x [L> 2D (fc) + L> 2 ^(fc)]^ | 2 (k) 



with 
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Here, we summed over N* central lattice sites. In the 
limit of an infinite lattice the maximum of moves 
towards k = with lim/^o U^)(k) — 1. Therefore, the 
total DDI potential for an infinite stack of BECs does 
not vanish anymore at k = (dashed line in the inset 
of Fig. 2). However, this is a pathological case because 
for any finite N s the total DDI potential vanishes at k = 
and our assumption of slowly varying wave functions 
breaks down towards the boundary. 



III. VALIDITY OF THE 2D MODEL 

In this section, we investigate the validity of the ef- 
fective 2D model for multilayered dipolar BECs intro- 
duced in Sec. II. To this end we computed ground 
states for the 3D GPE [Eq. (2)] [48], the coupled 2D 
model [Eq. (4)] [49], and the single mode 2D model 



[Eq. (12)] [50] using the normalized gradient flow (imag- 
inary time) method. For the time discretization we used 
backward Euler finite difference [50]. For the spatial dis- 
cretization we employed the sine pseudospectral [48] and 
the Fourier pseudospectral methods [40] for the 3D GPE 
and the 2D models, respectively. For the 3D compu- 
tation we assumed that the wave function vanishes at 
the boundaries. We integrated the 3D ground states 
over the individual lattice sites to find the N s densities 
IV>f(p)| 2 = Isle-i/2) dz\iP(r)\ 2 . To determine the valid- 
ity of the 2D model we compared the 2D ground states 
i/je(p) to -0| D (p). Using the single mode approximation 
reduced the computation times drastically: typically to 
less than a minute, compared to 2-3 hours for the cou- 
pled equations and ~ 1 day for the 3D GPE. In this 
section we only consider polarization in the x-z plane, 
that is d = (sin$,0,cos$) (cf. Fig. 1). Because the ex- 
ternal potential is radially symmetric, this simplification 
corresponds to choosing the transverse projection of the 
polarization direction as the x axis. 

To compare the axial profiles of the coupled 2D and 
3D ground states we computed the relative particle num- 
bers in each lattice site. Because of the long range of 
the DDI, we observe fairly pronounced boundary effects 
in the 3D computations for strong dipolar interactions 
g d ~ g. For this reason we omit the Nt outermost lat- 
tice sites in the overall normalization. Then the relative 
number of particles in site £ for the 2D model is given by 

N e = Jd 2 p\Mp)\ 2 /E^: +Nb Jd 2 p\Mp)\ 2 (the rel- 
ative particle number 2V? D for the 3D GPE follows by 
replacing |V^| 2 with |"0| D | 2 )- Figure 3 shows the parti- 
cle number difference (Nf D — N^)/Nq D relative to the 
particle number at the central lattice site. Although the 
number difference varies slightly over the central lattice 
sites, the difference between the GPE and the 2D model 
Eq. (4) remains smaller than 4% and 1% for the two pa- 
rameter sets, respectively. 

Next we compared the density profiles of the central 
lattice site ^(p)! 2 for the coupled and single mode mod- 
els with \ifjQ D (p)\ 2 ■ The sum of the densities of the cou- 
pled 2D and the total density of the 3D GPE are nor- 
malized to a function proportional of the total particle 
number Af(N). However, in the single mode approxima- 
tion we only consider a single wave function which has, 
consequently, a normalization less than Af. If the BEC 
density were the same in all layers, the normalization of 
this single wave function would be N /\/W s . Because the 
density varies slightly across layers, instead we chose to 
normalize the single mode density to the particle number 
in the central layer of the GPE. The ground state densi- 
ties for various DDI strengths and polarization angles are 
shown in Fig. 4. We find that both the coupled and single 
mode models describe the ground state well for any polar- 
ization. We only observe a slight difference between the 
models for strong DDI on the order of the contact inter- 
action and parallel polarization (top left panel in Fig. 4). 
This means that even the single mode approximation de- 
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FIG. 3. (Color online) Relative particle number difference 
between GPE ground state and the 2D model [Eq. (4)] for 
individual lattice sites. The particle numbers are relative to 
the particle number in the central layer Nq D (bars). The 
discs indicate the particle number difference in the 2D model 
relative to the central site (right axis label). The parameters 
are N s = 61 lattice sites with Vb = 20E r , E r /huj = 60, and 
g = 100V2£V//kj7r 2 . 



™ 



-2 
2 



-2 
2 



2A — i£ 

9 20 

♦ 



3d — I 
* 



1 ' 1 1 1 


1 1 1 1 1 


1.1,1 


1,1,1 


1 ' 1 1 1 




1 1 1 1 1 



1 1 1 1 1 



1 1 1 1 1 

i.i.i 



.1? = 



■i? = 



-2 2 -2 2 

x (units of «o) x (units of ao) 



I 



0.025 0.125 

\ipo\ 2 , |^o D | 2 (arbit. units) 



0.225 



scribes the ground state of the multilayer dipolar BEC 
well. Its accuracy diminishes for strong DDI because the 
true densities vary sufficiently strongly over the central 
lattice sites. 



IV. INTERLAYER-DDI-INDUCED CHANGE OF 
THE ASPECT RATIO 



FIG. 4. (Color online) Ground state densities of the central 
lattice site for various DDI strengths and polarization angles. 
The filled surfaces are the projection of the central site of the 
GPE results, whereas the solid (dashed) contour lines are the 
ground states of the coupled (single mode) 2D equation (4). 
The plotted densities are all normalized to 1. The coupled 
and single mode results are almost indistinguishable except 
in the top left panel. The parameters are as in Fig. 3. The 
plots use the magnetic length ao = y^h/muj as length unit. 



The interlayer DDI can cause observable effects in mul- 
tilayered dipolar BECs. This becomes apparent from 
Fig. 2. The strength of the interlayer DDI is compara- 
ble to the strength of the intralayer DDI at wavelengths 
larger than 5. We expect that the anisotropy of the DDI 
for d > leads to a change in the aspect ratio of a quasi- 
2D dipolar BEC in the central layer of a stack of dipolar 
BECs. In this section, we investigate these effects nu- 
merically using the single mode approximation for the 
central layer. 

To determine the mean radii of the central layer first we 
computed ground state densities for a varying number of 
lattice sites at a constant normalization. We calculated 
the mean radii as 

R% = [ d 2 pa 2 \Mp)\ 2 , (a = x,y). (15) 



The aspect ratio of the central layer is then given by 
Ry/Rx- Magnetostriction causes the dipolar BEC to ex- 
pand along the polarization direction [3, 40]. Figure 5 
shows the aspect ratio as well as the individual mean 
radii of the BEC as a function of the number of lattice 
sites N s . The case N s = 1 corresponds to a single layer 
dipolar BEC. We observe that the interlayer DDI causes 
an additional reduction in the aspect ratio depending on 
the number of lattice sites and polarization angle. For 
perpendicular polarization the aspect ratio remains un- 
changed because the DDI is isotropic. However, the in- 
dividual radii decrease. We have also computed aspect 
ratios for a stronger lattice with Vb = 40E r and observed 
a similar dependence of the mean radii on N s . For this 
stronger lattice and $ = ir/4 the aspect ratio was closer 
to 1 and its change slightly smaller than at Vq — 20E r . 



6 



13 
3 



1.0 



0.8 



1.00 

0.95 
1.2 



1.0 



' — # = tt/2 


.,--»-- »' *" 
< • - 


». * 

i 


— ♦ « - 


. _ — . — . — _ 


i 

— ■» » «.. 


d = tt/4 





— *■—-—♦« 


i 


# = 












Ry/Rx 

Rii 



Ry 

Ry/Rx 

Rx, Ry 
Ry/Rx 



11 

Total number of lattice sites N s 



21 



ary state 4>e{t), that is, ipi(p,t) = e '^[-07 + £(p, t)]. 
We expand the perturbation in a plane wave basis as 
&(p,t) = (1/2tt) / ^q^e 1 ^^-^*) + v*^-^-^*)) 
and insert ipe(p,t) into Eq. (4). Here, w g are the exci- 
tation frequencies of quasimomentum q and u q ^, Vq£ are 
the mode functions in layer I. Keeping terms linear in 
the excitations and Vqi we find the Bogoliubov-de 
Gennes equations for perpendicular polarization 



WqU q £ = — U^i + v(g + 2g d )(Uql + Vql) 



(17) 



FIG. 5. (Color online) Mean radii and aspect ratio of the 
central BEC layer as a function of the number of lattice sites. 
The different panels correspond to different polarization an- 
gles. The interlayer DDI has a noticeable effect over sev- 
eral lattice sites. The lines are marked at the right and are 
only to guide the eye. The parameters are as in Fig. 3 with 
g d /g = 19/20. 



For perpendicular polarization the mean radii and aspect 
ratio were nearly indistinguishable from the top panel 
in Fig. 5. The DDI-induced change of aspect ratio has 
been observed in a single layer 52 Cr via time of flight 
expansion[f5, 51]. We suggest that the dependence of 
the aspect ratio on N s could also be observed via time 
of flight expansion. To observe the central layers, in this 
experiment the outer layers would have to be removed 
on a time scale short enough to suppress equilibration, 
e.g., with additional lasers focused on the outer layers. 
This is followed immediately by time of flight expansion 
of the BEC. The observable effect is largest for parallel 
polarization $ = n/2. 



V. BOGOLIUBOV EXCITATIONS 

In this section we investigate the influence of interlayer 
DDI on the excitation spectrum of a layered quasi-2D 
dipolar BEC. In particular, we consider local density fluc- 
tuations of the layered BEC and derive their Bogoliubov 
energy. Their Bogoliubov energy can assume imaginary 
values for suitable parameters, which indicates the on- 
set of a dynamical instability that leads to exponential 
growth of excitations. 

To determine the Bogoliubov energy we consider small 
perturbations around the ground state of Eq. (4). For 
simplicity we assume a vanishing transverse harmonic 
potential Vho = and homogeneous density v in each 
layer. For an optical lattice with N s sites v = 1/N S . A 
stationary state of the effective 2D GPE (4) is given by 
tpe(p,t) = tpe(t) = e _1Mt y / z/ with the chemical potential 

ti = [g-g d 0--3d 2 z )]v. (16) 

Now we add a local perturbation £e(p,t) to the station- 



JqVqe = — w q £ + v(g + 2g d )(vqi + Uqt) 



(18) 



Excitations in layer I are coupled to excitations in all lay- 
ers through the interlayer DDI. However, the interlayer 
DDI drops exponentially with the distance [cf. Fig. 2 and 
Eq. (9)]. Therefore, first we only take into account near- 
est neighbor interactions \£ — j\ < 1. Then the matrix of 
the system of Eqs. (17)— (18) becomes tridiagonal and can 
be solved for its eigenenergies. The resulting Bogoliubov 
energy Eb (q) = w q is determined by 



E%(q) 



2 



2(g + 2g d )v 



3gdvU2u(q) - l2gdvU 2 t> ' (<?) 



(19) 



Because U^(q) vanishes for zero quasimomentum, the 
speed of sound c = Iim g _>o 9Eb (q) jdq = \J~gv + 2g d v is 
not influenced by the interlayer DDI. Only the intralayer 
DDI increases the speed of sound via its zero momentum 
mode. 

Now we generalize the Bogoliubov energy in multilayer 
dipolar BECs to arbitrary polarization. After inserting 
the expansion of the 2D wave functions into Eq. (4) we 
find the squared Bogoliubov energy 



£l(q) 



2 



2[g - g d (l - M 2 z )]v 



+ 6&ii/W&(q) - 12«fei/| 



(20) 



Here, 



[(d x cos ip + d y sin ipY 



dlP&iq) 



m 

polar coordinates q = q(cos<p,simp). In general, this 
excitation energy is anisotropic but mirror symmetric 
around the polarization direction projected onto the x—y 
plane. Interestingly, the interlayer interaction always re- 
duces the Bogoliubov energy compared to the Bogoliubov 
energy of a dipolar BEC with only intralayer DDI. This 
means that interlayer DDI drives the BEC closer towards 
an instability regardless of the polarization direction. 
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We gain qualitative insight into instabilities by look- 
ing at the dipole-dominated regime with g/ gd 0. Set- 
ting g — in Eq. (20) we see that the contact inter- 
action terms becomes attractive for polarization angles 
d 2 z = cos 2 < 1/3. Because the DDI terms (last line) 
in Eq. (20) vanish at q = 0, this leads to imaginary Bo- 
goliubov energies at low quasimomenta q. The dipole- 
dominated quasi-2D BEC ground state is not stable in 
this regime. However, a repulsive s-wave interaction 
g > g<i(l — 3<i 2 ) prevents this type of instability. For 
repulsive contact interaction (cos 2 $ > 1/3) another in- 
stability of the dipole-dominated quasi-2D BEC occurs 
at nonzero quasimomenta. For perpendicular polariza- 
tion a sufficiently large negative intralayer DDI term in 
Eq. (19) (large g^y) compensates the positive free energy 
and local terms (first line). This leads to an instability 
in the cross-over regime from quasi-2D to 3D [9]. The 
interlayer DDI term in Eq. (19) shifts the instability re- 
gion to smaller quasimomenta. Because in the present 
article we only consider quasi-2D BECs, we refer to an 
upcoming article investigating instabilities in the 2D-3D 
cross-over regime [52]. 

Figure 6 shows the Bogoliubov energy Eq. (20) of a 
dipole-dominated quasi-2D multilayer BEC for three po- 
larization directions. For nonperpendicular polarizations 
the Bogoliubov energy becomes anisotropic with higher 
energies along the projected polarization direction. In 
Fig. 6 we observe the instability at low momenta for 
d = tt/2. The cuts in Fig. 6(b) show the development of a 
roton minimum at moderately large DDI strength. The 
interlayer DDI advances the development of this mini- 
mum to smaller values of g& compared to a single layer 
quasi-2D dipolar BEC. For comparison we also plot the 
Bogoliubov energy for 52 Cr in Fig. 6(b), where we as- 
sumed that the contact interaction has been reduced to 
g = via a Feshbach resonance [16]. The interlayer DDI 
strength of 52 Cr is too weak to influence the dispersion 
significantly. 



VI. CONCLUSION 

We showed that interlayer DDI in a multilayer stack 
of dipolar BECs markedly reduces the aspect ratio of the 
quasi-2D BEC in the central layer. The greatest change 
in aspect ratio occurs for parallel polarization. We sug- 
gested that this effect of the interlayer DDI is observable 
in time of flight image of the central layer. 

To simplify numerical computations we presented a 2D 
model for a stack of quasi-2D dipolar BECs created by 
a strong ID optical lattice and transversely trapped in a 
harmonic potential. Our model is based on a dimension 
reduction of the GPE assuming a Gaussian axial den- 
sity profile of the wave function in the individual layers. 
We derived effective intra- and interlayer DDI potentials 
for the resulting coupled quasi-2D BECs. For weak in- 
terlayer DDI we observed only small variations in the 
particle numbers per lattice site, which allowed us to de- 
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FIG. 6. (Color online) Bogoliubov energies for different polar- 
izations and DDI strengths. The polar plots in (a) are marked 
with the magnitude and angle of q. White areas mark un- 
stable regions, (b): Cuts through Bogoliubov energies at the 
polar angles indicated in (a). Solid lines include intra- and in- 
terlayer DDI, whereas dashed lines only include the intralayer 
DDI. The green line represents 52 Cr. The interlayer DDI does 
not influence high energies where the in-plane excitations be- 
come particle-like. Parameters are as in Fig. 3 with g = and 
f = l/10. 



rive a single mode approximation for the quasi-2D BECs 
in the central sites. This approximation reduces the nu- 
merical computation of mean-field ground states of this 
system from ~ 1 day to several seconds. The resulting 
ground states match the reduced ground states of the 3D 
GPE excellently up to moderately large DDI strengths. 
For large DDI strengths g e i — g we still found very good 
agreement at all polarizations. 

Finally, the interlayer DDI reduces the squared Bogoli- 
ubov energy, which influences the development of a roton 



minimum and possibly leads to an instability (imaginary 
energy) for large density or DDI strength. The excitation 
spectrum of local perturbations becomes anisotropic for 
nonperpendicular polarization. 
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Appendix A: Derivation of the effective 2D model 

In this appendix we present the derivation of the ef- 
fective 2D model for multilayered dipolar BECs in a 
ID optical lattice [Eq. (5)]. First we use the identity 
U dd (r) = -c dd [5(r)/3 + a dd (l/47r|r|)] to split the DDI 
into a local and nonlocal part [48, 53]. Then we in- 
sert ip(r,t) = e -1 */ 27 J2jipj(p,t)wj(z) with Wj(z) = 
(l/7r7 2 ) 1 / 4 e _ ^ 2_z ^ 2 / 272 into Eq. (2), where we approx- 



imate V (z) ~ 2r*Y,j( z ~ z ]) 2 - We multiply by wi(z) 
and integrate the resulting equation over z. Setting 
J dzwi(z)wj(z) = for £ ^ j and using J dzw 2 (z) = 1, 
/ dzwj(z) = l/^rhry 2 we find 



id t ipe = 



-oVi + Vh +g(l-e d d)|^| 5 



rpt + Vt. (Al) 



d yy and 



Here, Vj_ = d x 

= ~i9d J dzd 3 r'w e (z)d dd U 3D {r - r') 



x tf(P' ''WpiP' \t)ip q (p,i)w j {z')w p (z')wg(z). 

(A2) 



The kernel in Eq. (A2) fulfills V 2 £/ 3D (r) = -6(r) so that 
for any / = /(r) 



d„(U SD *f) = -f-V 2 x (U SD *f), 



(A3) 



where * denotes a convolution. We expand the direc- 
tional derivative in Eq. (A2) as ddd = dd ± d± + d 2 z d zz + 
2d z dd ±z with dj^ = (d x ,d y ). Applying Eq. (A3) to the 
convolution in Eq. (A2) yields 



^2 J dzd 3 r'i) q (p,t)wq(z)we(z) 



x {d d±d± - d 2 z V 2 ± + 2d z d d±z )U 3D (r - r') 
x ip](p' ',t)tp p (p' ',t)wj(z')w p (z') . 



(A4) 



The first term in Eq. (A4) contributes to the contact 
interaction, whereas the second term forms the nonlocal 
potential. 

The even kernel U^ en [Eq. (6)] is determined by the 
terms in Eq. (A4) with only radial derivatives. After 
inserting C/ 3 d and the Gaussians Wj into Eq. (A4), we 
need to solve the integral 



dzdz' 



e -[(z'-z : ,) 2 + (z'~z p ) 2 + (z-z q ) 2 + (z-z e ) 2 ]/2-y 2 

Air 2 j 2 y/(x - x') 2 + (y - y') 2 + (z - z') 2 ' 



(A5) 

We substitute C, = z — z' — {z q + zt — Zj — z p )/2, £' = 
z + z' — (z q + zg + Zj + z p )/2 in Eq. (A5) and integrate 
over The solution defines the even ker nel of the DDI 
potential with p = yj(x — x') 2 + {y — y') 2 



2(2tt) 3 /2 7 



d(- 



,-CW p -(*i P +^W 



2 



(A6) 

In Fourier space with k = k(coscp, sin<^) the deriva- 
tives dd ± d± — d 2 \/\ in Eq. (A4) become — k 2 [(d x cos <p + 
d y sirup) 2 — d 2 ]. We use the convention /(k) = 
(l/27r) j d 2 pf(p)e~ lk p for the 2D Fourier transform. 
With this normalization the convolution theorem is F[f* 
h] = 2nJ : [f]J-[h]. For radially symmetric f(p) — f(p): 
f(k) — J dppf(p)J (kp) with J the Bessel function. Us- 
ing this formula for the Fourier transform of the Eq. (A6) 
and multiplying by 27rfc 2 from the convolution and the 
Fourier transform of the derivatives in Eq. (A4) we find 



k 



x e 



/ 7 2 fc + (S qj + 6 tp )/2 
( 1 2 k-{5 qj +8 lp )/2 \ 

2S 2 p +2S 2 ql , + (S qj +S,_ p ) 2 
8 7 2 



(A7) 



For j = p = q = £ Eq. (A7) reduces to the intralayer 
DDI U2T>(k). Because of the exponential prefactor, terms 
where all j,p, q, £ are mutually unequal are strongly sup- 
pressed. Similarly, terms with q = j, p = £ and j ^ £ are 
exponentially suppressed. The remaining terms q = I, 
p = j, and j 7^ £ form the interlayer DDI kernel L^ven 
[Eq. (6)]. 

The odd kernel U 3 oAA [Eq. (7)] is determined by the 
term in Eq. (A4) with derivative dd ±z - Using d z {UsD * 
<?) = (dzll^u) * 9 we insert the derivative dzU^D into 
Eq. (A4). Then we need to solve the integral 



dzdz 



, (z - z')c-[( 2 '- z >) 2+ ( z '-^ 2+ ^- z «) 2+ ( z ^ 2 ]/ 272 

47T 2 7 2 [(a; - x') 2 + (y - y') 2 + (z - z') 2 } 3 / 2 

(A8) 

Following the steps for the even kernel we obtain the odd 
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kernel 



1 



dt C 



2(2tt) 3 / 2 7 7 V' ' 1 

e -C 2 /27 2 e -(^p+^)/47 2 



(A9) 



p2 + ( c + ^+vy 



3/2 ' 



The Fourier transform of J7^| [Eq. (A9)] multiplied by 
2irk from from the Fourier transforms of the convolution 
and the remaining radial derivative is given by 



/ 7 2 fc + (S qj + 5 ep )/2 
/ Vfc-(<^-+fr p )/2 \ 

2^ p +2^ + (^+^ P ) 2 
8 7 2 



(A10) 



Only terms with q — £, p — j are not exponentially sup- 
pressed in Eq. (A10). Hence, we recover U^ dd [Eq. (7)]. 

By combining Eqs. (A7) and (A10) with Eq. (A4) and 
neglecting the suppressed terms in the sum we recover 
the DDI potential Eq. (5) in Fourier space. 

For completeness we present an approximation of the 
spatial potential for multilayer DDI with arbitrary polar- 



ization direction. To obtain this approximation we take 
the limit 7 —> in Eqs. (A6) and (A9) treat the Gaussians 
in £ as approximations for the Dirac delta distribution: 



7^0 4-7T 



1 



1 Si 



1/2' 



(All) 
(A12) 



Again we neglect the exponentially suppressed terms. In- 
serting these kernels into Eq. (A4) and calculating the 
remaining derivatives we find 

Vf D (p) = 3 ^E/ d P' U 2u(P P'M(p',t)\ 2 (A13) 



with 
Ug{p) 



1 



5/2 



[P 2 + (1 - 3dM 



4^(W<^y (A14) 
-64%d x -p-3|d x -p| 2 ]. 

For the intralayer part j = I this approximation remains 
valid for p 7. We note that Eq. (A14) corresponds to 
the dimensionless DDI potential Eq. (1) projected onto 
2D planes separated by Sgj. 
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